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An effective field theory exists describing a very large class of biophysically interesting Coulomb 
gas systems: the lowest order (mean-field) version of this theory takes the form of a general- 
ized Poisson-Boltzmann theory. Interaction terms depend on details (finite-size effects, multipole 
properties, etc). Convergence of the loop expansion holds only if mutual interactions of mobile 
charges are small compared to their interaction with the fixed-charge environment, which is fre- 
quently not the case. Problems with the strongly-coupled effective theory can be circumvented 
with an alternative local lattice formulation, with real positive action. In realistic situations, with 
variable dielectric, a determinant of the Poisson operator must be inserted to generate correct 
electrostatics. Methods adopted from unquenched lattice QCD do this very efficiently. 
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1. Introduction 

A central paradigmatic problem of modern biophysics centers on the thermodynamic proper- 
ties of generalized Coulomb gases, by which is meant a system of both mobile and fixed electric 
charges (or more generally, particles with higher multipole moments as well), in a background (typ- 
ically water, but also electrically neutral interiors of protein polypeptides, in which fixed charges 
may be embedded) of variable dielectric. Direct molecular dynamics simulations of Coulomb gas 
problems are impractical for large systems, as the long-range character of the Coulomb interaction 
means that the electrostatic energy of every pair of particles in the system has to be computed. In 
systems with uniform dielectric, this problem can be ameliorated by Fourier (Ewald) techniques, 
but in the general case, the dielectric "constant" in fact varies spatially (and also dynamically, in 
the course of the simulation, if macroions or polymer components are allowed to change confor- 
mation), making Fourier methods impractical. 

Of course the screening effects in a homogeneous medium with mobile charges have been 
understood for a long time on the basis of Debye-Hiickel theory |[T]], provided the charge concen- 
trations are not too high. Here one starts from a Poisson-Boltzmann equation which summarizes 
the mean-field effects of the mobile ions on each other (and with any fixed charges). Starting in the 
1940s, a great deal was accomplished with a linearized version of this equation (the DLVO method 
introduced by Derjaguin, Landau, Verwey and Overbeek |Q]). However, such methods do not pro- 
vide any intrinsic procedure for systematic improvement of the mean-field result, and frequently 
fail completely in the regime of high concentrations. 

A more general formalism for dealing with Coulomb gases under very general circumstances 
was introduced in the early 90's the initial emphasis was in dealing with systems of fixed 
charged macroions surrounded by a gas of small mobile ions. The grand canonical partition func- 
tion for such systems can be converted into a path-integral regularized on a spatial lattice, and the 
saddle-point expansion of the functional integral then leads to a (discretized) Poisson-Boltzmann 
equation which can be rapidly solved numerically. Moreover, the basic technique allows for 
straightforward generalization to systems with short-range repulsive forces between the mobile 
ions 10], multipolar ions [^, charged polymer interactions [||, 0, |]|, among others. In all of these 
cases the higher order fluctuation ("loop") effects are clearly defined, and the leading corrections 
to the mean-field result computable. Despite the formally attractive nature of this effective field 
approach, there still remains the difficulty that in many cases these fluctuation corrections are very 
large, so a perturbative saddle-point expansion does not yield useful results. Furthermore, the effec- 
tive action for these theories is always complex, so direct Monte Carlo simulation of the functional 
integral is impractical due to the infamous sign effect. 

Recently, Maggs and collaborators suggested [^] an alternative approach to Coulomb gases in 
which the long-range Coulomb interaction is localized by writing the path integral for the parti- 
tion function in terms of local electric field variables (essentially, one takes the Hamiltonian path 
integral for finite temperature Maxwell electrodynamics and neglects magnetic terms). A consid- 
erable amount of work has now been devoted to streamlining and improving the efficiency of this 
approach [^, 11]. The original algorithm of Maggs et al. only handles systems of uniform dielec- 



tric, however, for which one may argue that Fourier accelerated molecular dynamics simulations 
are competitive. For systems in which the dielectric medium is dynamical, the Maggs et al. func- 



2 



Lattice Field Theory Methods in Modem Biophysics 



Anthony Duncan 



tional integral produces a spurious interaction force between the particles, which must be removed 



to obtain the correct electrostatic energy. Recent work \ \12[ |13| ] has shown that this can be done 
exactly by introducing the determinant of the generalized Poisson operator into the path integral, 
in complete analogy to the way the determinant of the quark Dirac operator must be introduced in 
the path integral of unquenched QCD. Numerical simulations can be performed using a variety of 
methods imported from lattice QCD: in the following, we describe results for the structure factor 
of a dielectric plasma obtained using the Liischer multiboson method. 

2. Functional Formalism for Coulomb Gases 

Although the Coulomb interaction is long-range, it has the special, and extremely useful, fea- 
ture that it is the Green's function for a local operator, the Laplacian. This makes it possible to 
remove the nonlocal Coulomb interaction by introducing a local auxiliary field- 

^-ljdrdr'p{r)p{r')/\r-r'\ 



Here p(r) is the charge density: to be specific, imagine a system with stationary charges 
(charge Qje, locations Rj) and mobile simple ions (charges ±e, locations ric,ri) so the charge 
density is 

p(r) = -e'£Qj5{r-Rj)+e'£5{r-r,) 

j k 

The full Hamiltonian in a typical Coulomb gas consists of a Coulomb energy term, plus single- 
particle exclusion potentials 



2e J Ir-fl 

+ lLyi7k) + ^ZV{ri) 
P k Pi 

It is most convenient to work in the grand-canonical framework, so we introduce chemical 
potentials and the partition function becomes 



Y[drkdrie 



In addition to localizing the Coulomb potential, the auxiliary field representation factorizes the 
partition sums: 

■ dfkdri 

• exp ( - £ (r/ ) + V (r/) ) + j3 (At+«,t + At-«/ ) ) 



£ /^exp(£(/a(r,)-V(r,))) 



: exp {el^^^ J dre^^^-^^+eP^"- J dre^-''^-^^] 



3 



Lattice Field Theory Methods in Modern Biophysics 



Anthony Duncan 



and we obtain an exact path integral representation for Zgc- 



Z, 



^Z(^)exp(- J x{r)Ax{r)dr) 

exp ( I ir+e'^-^ + Y^e-'x-^)dr) 
exp{-iY,QjX{Rj)) 



J 

&Xir)eW-Six) (2.1) 

where a = j^^^ , 7± = e^^^^ . Although we have written the path integral in terms of fields defined 
for continuous space, one needs to regulate the preceding (and following) formulas on a spatial 
lattice to have a well-defined theory. 

2.1 Mean-field Theory and Loop Expansion 

The grand-canonical partition function Zgc = / &x^'^V~^{.X) has a complex saddle-point of 
the Hubbard-Stratonovich type: 

X = -4 + ^ 



5x{r) 

The saddle-point condition is the Poisson-Boltzmann equation for the system: 

j 

In general, this is solved rapidly and accurately by relaxation techniques. Evaluating the action at 
the saddle-point gives the Thermodynamic Potential S = log Zgc in leading order: 



a [ - - 



- I^QjHRj) (2.2) 

j 

In practice, all spatial integrals and differential operators are regulated on a spatial lattice (as men- 
tioned previously). 

2.2 Systematic Loop Expansion of Thermodynamic Potential S 



The mean-field approximation ( |2.2[ ) is only useful if the fluctuation corrections are small: in 
fact, if they are, we have a systematic procedure for improving the mean-field result. Expanding 
around the saddle-point (mean-field) solution, x = —i^ + ^ '■ 



— -^mf ~l~ -^loops 
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where 

Sloops = log(| ^^exp-(5o(^)+5int(^))) 

The mean-field action is 

and the fluctuation terms (through order E,^) are 

The loop (i.e. fluctuation) corrections involve a propagator defined as 

D{r,r) =< r\ {-aA + y+e^-^ + 7_e-^-^)-' |r' > 

Unlike the situation in conventional field-theoretical perturbative expansions, in typical Coulomb 
gas applications this propagator cannot be written down analytically, as the system involves fixed 
charges in essentially arbitrary locations (so Fourier methods fail). Instead, the propagator is com- 
puted (on a lattice) numerically. Once this is done, 1-loop, 2-loop, etc. corrections to the thermo- 
dynamic potential follow immediately: 

S(i) = — !-logdetD 
2 

The formalism described above has been applied to numerous problems, including the thermo- 
dynamics of charged polymers in electrolytes [^] and the partitioning of charged polymers between 
cavities 

2.3 Extensions of the Formalism 

The path integral formalism for Coulomb gases outlined above is extremely flexible: one may 
easily modify it to include 

1. Non-Coulomb pairwise interactions (such as repulsive Yukawa) of the mobile charges [Q]. 

2. Higher multipoles (dipole, quadrupole, etc) on the mobile charges 

In order to introduce a repulsive short-range core interaction between the mobile charges (in 
addition to the long range Coulomb piece) we consider |Q] a Coulomb/Yukawa Gas, with the inter- 
action energy of mobile charges given by: 
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+ ^ / drdr'p,{r) p,{r') 

^ J \r — r\ 

j k 

I 

Ps(r) = £5(r-r;t)+£5(r-r/) 

k I 

One then introduces a second auxiliary field for the Yukawa component: 

exp{— ^ / drd?ps{r f ' ^ ' p.(P)} 

= j ^'^('')^^P{^^ / ¥i^-^^)¥dr 
+ iVc I W)p,{r)dr} 



Integrating over mobile charge positions as before yields an equivalent theory in terms of two 
auxiliary fields i/a: 



2gc = j &X{r) j ^V^(r)exp{g|^ j x^xdr 



with 



Smi{X,¥) = j dre 



-i:JeX{r)+'\^V{r)-V(f) 



e 



11 

1 

Similarly one may derive an extension of the formalism to deal with a multipolar gas, in 
which mobile charges also characterized by higher multipoles (dipole, quadrupole, etc): 

j 
j 
j 



6 



Lattice Field Theory Methods in Modem Biophysics 



Anthony Duncan 



The total electrostatic energy for such a gas isgiven in terms of the effective charge density: 

Peff(r) ^ p(r) - daP"{r) + badbQ"\7) + .. 

o 



■A' 

As usual, we can localize the long-range interactions with an auxihary field: 

-B f ^ 1 
exp{— / i/^r^/Vpeff(r)— — peff(r')} 

^Dxir)exp{ J i:^X^X + iXir)pis{r))dr} 

As an example, consider a gas of mobile permanent dipoles, \pj\ = p (Q"^, .. = 0) The config- 
uration integral for a single dipole involves the average over orientations: 

leading to derivative interactions in the effective field action: 

SmtU) = el^^^ J e^'^^^'^-^^^Ur 

J p\yx{'-)\ 

so we are led to a modified Poisson-Boltzmann equation, in which the mobile dipoles provide an 
effective spatially varying dielectric, together with a correspondingly modified loop expansion. 

3. Weakly/Strongly Fluctuating Coulomb Gases 

In a renormalizable local field theory like QCD, there is a natural dimensionless coupling 
(typically, the running coupling at momentum scales relevant to the process under consideration) 
which provides an expansion parameter for the saddle point expansion corresponding to covariant 
perturbation theory. In the case of the effective field theories discussed above for Coulomb gas 
problems, the validity of a perturbative loop expansion around the mean-field (Poisson-Boltzmann) 
theory depends on the ratio of two length scales: 

1. The Bjerrum length /g = distance between two mobile charges such that pair electrostatic 
energy = ^gT. 

2. The Gouy-Chapman length fi, which depends in a more complicated way on the geometry 
of the fixed charges relative to the mobile ones. 
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Figure 1: Coulomb gas near charged plate 



As a simple example (see [14]), consider a gas of mobile ions of charge —qe near a charged 
plate of charge density +oe (see Figure 1), where the whole system is electrically neutral. For this 
problem, the Bjerrum and Gouy-Chapman lengths are 



h 



In this case, the Gouy-Chapman length /i corresponds to the distance from the plate at which an 
isolated mobile ion has electrostatic energy ksT. The ratio E plays the role of the perturbative 
expansion parameter for this system, as can be seen by rewriting the Hamiltonian for the system 



1 



-H 



Ib 



j<k Ki ■ 



"7 + L^i 

rk\ q t 



in terms of distances rescaled to the Gouy-Chapman length, rj = pLry. 



ksT 



Note that the parameter E oc is extremely sensitive to the valence q- it often happens that we go 
from weak to strong fluctuations when q goes from 1 to 2 (monovalent to divalent ions) ! 

Unfortunately, in many interesting cases E > > 1 and the perturbative loop expansion breaks 



down. Direct Monte Carlo simulation of the path integral Eq.([2.1[) is not feasible:the action is com- 
plex and phase oscillations result in unmanageably large fluctuations (the infamous sign problem!). 
For such strongly coupled Coulomb gases (as for intrinsically nonperturbative field theories such as 
QCD in the infrared), we must resort to numerical simulation techniques, but clearly one needs an 
alternative formulation where the effective action is real, allowing the application of conventional 
Monte Carlo techniques. 
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4. Local Gauge Theory Approach to Coulomb Gases 

Start from the Hamiltonian path integral for abehan gauge theory with external point charge 
sources: however c oo and magnetic effects and time-dependence are ignored (as well as spatial 
variation of the dielectric). Then the canonical partition function takes the form: 

f ^ ^ ^ ^ 1 

Z= mdrim{r)lld{V-E--p{r)) 

J i=l r ^ 

i I 

The transverse, or curl part of the electric field variable decouples from the charged particle dy- 
namics via the Helmholtz decomposition: 

E = V(j) + V xA 
J drP = j Jr|V0|2 + j dr\VxA\^ 

The unphysical curl part of the electric field decouples from the gradient part: only the latter sees 
the charge density p , so that the charges couple to yield the correct electrostatic energy. In practice 
the functional integral is regulated on a spatial lattice, with mobile charged particles associated 
with sites and electric field values with links (see Figure 2). 
Rescaling to lattice units, the Hamiltonian becomes: 

and the Gauss' Law constraint takes the form: 

for the sum of outgoing link fields from any site containing a charged particle of charge qie. 
A simulation algorithm for this system is easily devised: 

1. Pick starting lattice locations (randomly) for the N particles of charge qi,i = l,..N. Then 
solve Gauss' Law for these fixed charge locations to obtain a starting configuration of electric 
link field variables satisfying the Gauss constraint. 

2. Update the electric fields by shifting all link variables along a complete set of independent 
closed paths by constant shifts. In the simplest version of this, one simply considers all 
plaquettes (unit squares) on the lattice, shifting the 4 link fields ordered around the plaquette 
by a random a,— (Vax < cc < (Xmax- This preserves Gauss' law. In fact, for systems with 
constant dielectric, one may use [ pTj ] fast Fourier transform methods to effect a global update 
of the electric fields, completely eliminating autocorrelations. 
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Figure 2: Electric Field Link Variables 

Update particle locations by visiting in turn every site n containing a charged particle of 
charge qj. A particle move to the neighboring site n + jX in a random direction /i is then 
considered, where the particle move is accompanied with a shift of the electric field Ei on 
the link I = {n + jl) 

Ei->Ei- qi 



to maintain Gauss' law. In practice [10], we need to couple particle moves to field changes 



on several neighboring links to get reasonable acceptance rates. 

A classic example of a strongly coupled Coulomb gas is the case of two like charged plates (see 
Figure 3), with mobile counterions (rendering the system overall neutral) between the plates. For 
appropriate choice of the charge density a on the plates, the system can be switched from weakly 
to strongly coupled simply by doubling the valence of the counterions (i.e. half as many doubly 
charged ions). For singly charged ions one finds, using the local gauge simulation techniques de- 
scribed above, a positive pressure, as shown in Figure 4 (also shown is the dependence on lattice 
discretization for the same physical size system). Mean field theory, along the lines discussed pre- 
viously in Section 2, always gives a repulsive plate interaction of this form. On the other hand, if 
we double the charge on each ion, leaving the plate charge density fixed (and halving the number of 
counterions), one then finds an attractive force between the like charged plates (negative pressure), 
implying a complete breakdown of the saddle-point expansion of Section 2. The result of a nu- 
merical simulation of this system using the local gauge method is shown in Fig. 5 (for details, see 



[10]), where the total pressure (solid line) is seen to be negative. These results agree with explicit 



(and painful) molecular dynamics simulations [ 15] of the like-charged-plate problem. 



5. Inhomogeneous Dielectric Effects in Coulomb Gas Problems 

In reahstic situations, the dielectric field e is NOT constant, but varies from 2-6 in the interior 
of macromolecules to 80 in the surrounding medium (water) 

Dir) = -e{r)y^{r)+y y.A{f)=b\\r)+rP{r) 
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Figure 3: Coulomb gas between like charged plates 
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Figure 4: Plate pressure for univalent ions 
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Figure 5: Plate pressure for divalent ions 
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Maxwell's equations imply = but if this constraint is not explicitly imposed: 

^2 r\\\(^\l r ntr/'^^^2 



:ir) J e(r) J e(r) 
Ignoring the irrotational constraint, if we procede as previously: 

Z' = / YldriYlm{r)6{V-D-p{r)) 

.g-f/drD7e(r) 

= / Y[driY[^D\\{r)m''{f)5iV ■ D\\ - p(r)) 

•' !=1 r 

. f /.fi^^ll (r)7e(r)g- f /cirS''(r)Ve(r) 

The integral over D"^ gives a spurious unphysical dependence on e{r)- which can change dynami- 
cally in the course of the simulation (e(r) depends imphcitly on rf). Regularizing the path integral 
with a spatial lattice, the exact form of the spurious e-dependence can be revealed simply by turning 
off all free charges so the contribution from the longitudinal part vanishes: 



^ w nil 



n nil 



n nfi 



Evidently, in order to eliminate the spurious term, we need to include in the path integral: 

^-i(e) = ]^e-5^"Mi°g(^M)det+5(-A^£^A^) 

Including a positive power of the determinant of a local operator in a path integral is precisely the 
problem we face in unquenched QCD ! 

5.1 Eliminating transverse contributions with Liischer multiboson fields 

In the Liischer approach to unquenched QCD, one begins with a polynomial approximation to 
1 /s in the interval [0, 1]; a convenient choice is: 

1 ^ 



s 



^i■k = :^(l + 5)(l-cos 



2' 2Af+l' 

^^ = ^^"^2ivTT 
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Figure 6: Structure Function for Neutral Dielectric Plasma 



implying for the determinant (scaling the spectrum to [0, 1]) a corresponding approximation: 

N 

n 

k=i 

This leads to the corrected form for the path-integral: 



det+5(^) ~ ndet-5((^-/ii)2 + vf^ 



Z= j dridD„^d<plt'^5{A^Dn^-pn 



The approximation only adequately describes the low eigenstates if one uses a sufficiently large 
number N of auxiliary fields, which limits its usefulness in QCD, where the density of low eigen- 
modes of the Dirac operator is large, due to spontaneous symmetry breaking. Fortunately, in bio- 
physical applications, the IR spectrum of the Poisson operator is sparse (unlike QCD) and we can 
get away with a small number of boson fields (as we shall se shortly, often N=A is adequate). 

A useful testbed for the multiboson implementation of a Coulomb gas with dynamical di- 
electric effects is the dielectric plasma, in which mobile particles (either neutral or charged) are 
assigned a dielectric constant different from that of the ambient medium. It is convenient to asso- 
ciate the dielectric field with links, each link being given a value depending on whether a particle 
is present at either, neither or both of the end sites of the link (for further details, see [|l^, 13]). An 
appropriate observable is the structure factor S(ci) = / dfe"i'^ < p(r)p(0) >: for neutral particles, 
this must be ^-independent. As a first test, consider a model with 1000 neutral dielectric particles, 
on a 16^ lattice, with a ratio of particle to background dielectric given by £part/£bg=0.2. 

As we see in Figure 6, once the multiboson fields are included, the measured structure factor 
is indeed flat, to a very good approximation. A somewhat more realistic model is the 1 component 
charged dielectric plasma, on a 32^ lattice, with 8000 particles, and £part/£bg=0.05. In Figure 7 
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Figure 7: Structure Function for Charged Dielectric Plasma 
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Figure 8: Dependence on number of multiboson fields 

we show the results for such a plasma, with and without the contribution from the multiboson 
fields. Evidently, the results are qualitatively incorrect without taking into account the induced 
effects in the transverse electric field. In Figure 8, we show the dependence on the number of 
multiboson fields used. As remarked previously, it is remarkable how few fields are needed to 
achieve reasonable accuracy in this model. 

6. Conclusions 

The last 15 years has seen an extensive development in the theory of the thermodynamics 
of Coulomb gases, vastly extending the scope of treatable systems beyond the simplest cases to 
which elementary Debye-Hiickel theory is applicable. The new methods treat a combined system 
of charged (and even multipolar) particles interacting with fields, and with both particles and fields 
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realized on a discrete spatial lattice. Formally, the functional formalism for the grand-canonical 
partition function offers the greatest generality, but is restricted in usefulness to weakly coupled 
Coulomb gases. The local gauge formulation introduced by Maggs and Rossetto [^, and recently 
generalized to systems with dynamical dielectric by Coalson, Duncan and Sedgewick [ |l2| ] allows 
the efficient simulation of strongly coupled Coulomb gases in a local formalism, circumventing the 
nonlocal Coulomb interaction which complicates, and frequently renders intractable, the molecular 
dynamics approach. 
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